REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources  gathering  and  maintaining  the - 

data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information  inducing  suggestions  for  reducing 
this  burden  to  Department  of  Defense.  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188).  1215  Jefferson  Davis  Highway  Suite  1204  Arlington  VA  22202-  * 

4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collodion  of  information  if  it  does  not  disdav  a  currentlv 
valid  OMB  control  number  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS.  ^  y  CU™nuy 

1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

29-Sep-2006  REPRINT 

3.  DATES  COVERED  (From  -  To) 

4.  TITLE  AND  SUBTITLE 

DETERMINATION  OF  CRUST  AND  UPPER  MANTLE  STRUCTURE  BENEATH 

AFRICA  USING  A  GLOBAL  OPTIMIZATION  BASED  WAVEFORM 

MODELING  TECHNIQUE 

5a.  CONTRACT  NUMBER 

FA87 1 8-04-C-00 1 4 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

62601F 

6.  AUTHOR(S) 

Jay  Pulliam,  Mrinal  K.  Sen,  and  Abhijit  Gangopadhyay 

5d.  PROJECT  NUMBER 

1010 

5e.  TASK  NUMBER 

SM 

5f.  WORK  UNIT  NUMBER 

A1 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Texas  at  Austin 

4412  Spicewood  Springs  Rd.,  Building  600 

Austin,  TX  78759-8500 

8.  PERFORMING  ORGANIZATION  REPORT 

NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory 

29  Randolph  Road 

Hanscom  AFB,  MA  01731-3010 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFRL/VSBYE 

11.  SPONSOR/MONITOR'S  REPORT 

NUMBER(S) 

AFRL-VS-HA-TR-2006- 1 160 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  Distribution  Unlimited. 

13.  SUPPLEMENTARY  NOTES  “  - - - 

Reprinted  from:  Proceedings  of  the  28th  Seismic  Research  Review  -  Ground-Based  Nuclear  Explosion  Monitoring  Technologies,  19-21 
September  2006,  Orlando,  FL,  Volume  I  pp  196  -  208. 

14.  ABSTRACT  During  the  initial  years  of  the  project,  a  waveform  modeling  program  based  on  the  reflectivity  method  and  incorporating  a  global  optimization  algorithm,  was  developed. 

Assuming  a  one-dimensional,  isotropic,  layered  Earth,  the  method  computes  synthetic  seismograms  by  distributing  independent  computations  for  different  ray  parameters  across  multiple  processors, 
providing  nearly  a  linear  decrease  in  computation  time  with  the  number  of  processors.  During  the  inversion  process,  the  method  implements  a  global  optimization  algorithm  using  very  fast  simulated 
annealing  (VFSA).  Global  optimization  allows  for  searching  within  a  broad  spectrum  of  models  in  order  to  find  the  global  minimum,  and  eliminates  dependency  on  the  starting  model.  In  particular, 
VFSA  is  advantageous  because  it  allows  for  larger  sampling  of  the  model  space  during  the  early  stages  of  the  inversion  process,  and  much  narrower  sampling  in  the  model  space  as  the  inversion 
converges  and  the  temperature  decreases,  while  still  allowing  the  search  to  escape  local  minima.  Additionally,  the  ability  of  VFSA  to  perform  different  perturbations  for  different  model  parameters 
allows  for  individual  control  of  each  parameter  and  the  incorporation  of  a  priori  information.  The  method  also  uses  statistical  tools  such  as  computation  of  the  posterior  probability  density  (PPD) 
function  and  the  parameter  correlation  matrix,  so  as  to  evaluate  the  uniqueness  of  a  particular  model  obtained  by  searching  the  model  space,  and  the  trade-offs  between  individual  model  parameters 
therein.  The  waveform  modeling  approach  has  now  been  applied  to  ground  truth  data  recorded  in  twelve  permanent  broadband  seismographic  stations  spanning  the  African  continent.  The  waveforms 
modeled  included  S,  Sp,  SsPmP,  and  shear-coupled  PL(SPL)  phases  from  deep  (200-800  km)  earthquakes  of  magnitude  5.5  to  7.0  located  at  distances  of  30  to  80  degrees.  Earthquake  data  that  satisfy 
the  search  criteria  described  above  are  sparse,  however.  Seventeen  events  meeting  the  above  criteria  were  available  for  analysis.  Forward  computation  of  the  synthetic  seismograms  was  performed 
using  models  from  independent  receiver  function  studies  (wherever  available)  as  an  initial  velocity  model  for  the  crust.  Preliminary  reference  earth  model  (PREM)  was  used  for  mantle  velocities 
below  the  receiver  function  models.  The  global  optimization  algorithm  with  several  hundreds  of  iterations  was  then  implemented.  This  algorithm  perturbs  the  crust  and  upper  mantle  velocity  structure 
(up  to  1 00  km  depth)  within  user-specified  bounds  to  generate  synthetic  seismograms  that  best  fit  the  data.  The  results  show  that  the  final  velocity  models  obtained  from  the  application  of  this 
method  for  various  stations  in  the  African  continent  correlate  well  with  those  obtained  from  receiver  function  techniques.  In  particular,  more  detailed  crust  and  upper  mantle  structure  was  obtained  for 
seismographic  stations  located  in  north  and  west  Africa,  where  prior  studies  are  sparse.  Broad  tectonic  observations  were  supported  by  the  models  generated  using  this  method;  specific  refinements  of 
the  models  for  some  stations  are  underway.  Overall,  the  method  is  successful  in  determining  the  crust  and  upper  mantle  structure  of  the  African  continent.  Notable  advantages  of  this  method  include 
ability  to  simultaneously  model  all  the  waveforms,  assess  uncertainties  in  model  parameters,  finding  a  range  of  acceptable  models  that  explain  the  data,  and  avoiding  a  dependence  of  the  final  model 
on  the  starting  model.  Furthermore,  it  differs  from  other  commonly  used  methods  in  the  fact  that  it  provides  a  direct  measurement  of  P  and  S  wave  velocities  simultaneously.  Also,  the  use  of  the  SPL 
phase  wherever  available  improves  constraints  on  the  models  of  the  lower  crust  and  upper  mantle. 

15.  SUBJECT  TERMS  - - 

Seismic  velocity.  Seismic  propagation 

a.  REPORT 

UNCLAS 


b.  ABSTRACT 

UNCLAS 


c.  THIS  PAGE 

UNCLAS 


OF  ABSTRACT 

SAR 


OF  PAGES 

13 


Robert  J.  Raistrick 


19b.  TELEPHONE  NUMBER  (include  area 
code) 

781-377-3726 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  239.18 


28th  Seismic  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


DETERMINATION  OF  CRUST  AND  UPPER  MANTLE  STRUCTURE  BENEATH  AFRICA  USING  A 
GLOBAL  OPTIMIZATION  BASED  WAVEFORM  MODELING  TECHNIQUE 


Jay  Pulliam,  Mrinal  K.  Sen,  and  Abhijit  Gangopadhyay 


University  of  Texas  at  Austin 
Sponsored  by  Air  Force  Research  Laboratory 
Contract  No.  FA87 1 8-04-C-00 1 4 


ABSTRACT 


During  the  initial  years  of  the  project,  a  waveform  modeling  program  based  on  the  reflectivity  method  and 
incorporating  a  global  optimization  algorithm,  was  developed.  Assuming  a  one-dimensional,  isotropic,  layered 
Earth,  the  method  computes  synthetic  seismograms  by  distributing  independent  computations  for  different  ray 
parameters  across  multiple  processors,  providing  nearly  a  linear  decrease  in  computation  time  with  the  number  of 
processors.  During  the  inversion  process,  the  method  implements  a  global  optimization  algorithm  using  very  fast 
simulated  annealing  (VFSA).  Global  optimization  allows  for  searching  within  a  broad  spectrum  of  models  in  order 
to  find  the  global  minimum,  and  eliminates  dependency  on  the  starting  model.  In  particular,  VFSA  is  advantageous 
because  it  allows  for  larger  sampling  of  the  model  space  during  the  early  stages  of  the  inversion  process,  and  much 
narrower  sampling  in  the  model  space  as  the  inversion  converges  and  the  temperature  decreases,  while  still  allowing 
the  search  to  escape  local  minima.  Additionally,  the  ability  of  VFSA  to  perform  different  perturbations  for  different 
model  parameters  allows  for  individual  control  of  each  parameter  and  the  incorporation  of  a  priori  information.  The 
method  also  uses  statistical  tools  such  as  computation  of  the  posterior  probability  density  (PPD)  function  and  the 
parameter  correlation  matrix,  so  as  to  evaluate  the  uniqueness  of  a  particular  model  obtained  by  searching  the  model 
space,  and  the  trade-offs  between  individual  model  parameters  therein. 

The  waveform  modeling  approach  has  now  been  applied  to  ground  truth  data  recorded  in  twelve  permanent 
broadband  seismographic  stations  spanning  the  African  continent.  The  waveforms  modeled  included  S,  Sp,  SsPmP, 
and  shear-coupled  PL  (SPL)  phases  from  deep  (200-800  km)  earthquakes  of  magnitude  5.5  to  7.0  located  at 
distances  of  30  to  80  degrees.  Earthquake  data  that  satisfy  the  search  criteria  described  above  are  sparse,  however. 
Seventeen  events  meeting  the  above  criteria  were  available  for  analysis.  Forward  computation  of  the  synthetic 
seismograms  was  performed  using  models  from  independent  receiver  function  studies  (wherever  available)  as  an 
initial  velocity  model  for  the  crust.  Preliminary  reference  earth  model  (PREM)  was  used  for  mantle  velocities  below 
the  receiver  function  models.  The  global  optimization  algorithm  with  several  hundreds  of  iterations  was  then 
implemented.  This  algorithm  perturbs  the  crust  and  upper  mantle  velocity  structure  (up  to  -100  km  depth)  within 
user-specified  bounds  to  generate  synthetic  seismograms  that  best  fit  the  data.  The  results  show  that  the  final 
velocity  models  obtained  from  the  application  of  this  method  for  various  stations  in  the  African  continent  correlate 
well  with  those  obtained  from  receiver  function  techniques.  In  particular,  more  detailed  crust  and  upper  mantle 
structure  was  obtained  for  seismographic  stations  located  in  north  and  west  Africa,  where  prior  studies  are  sparse. 
Broad  tectonic  observations  were  supported  by  the  models  generated  using  this  method;  specific  refinements  of  the 
models  for  some  stations  are  underway. 

Overall,  the  method  is  successful  in  determining  the  crust  and  upper  mantle  structure  of  the  African  continent. 
Notable  advantages  of  this  method  include  ability  to  simultaneously  model  all  the  waveforms,  assess  uncertainties  in 
model  parameters,  finding  a  range  of  acceptable  models  that  explain  the  data,  and  avoiding  a  dependence  of  the  final 
model  on  the  starting  model.  Furthermore,  it  differs  from  other  commonly  used  methods  in  the  fact  that  it  provides  a 
direct  measurement  ot  P  and  S  wave  velocities  simultaneously.  Also,  the  use  of  the  SPL  phase  wherever  available 
improves  constraints  on  the  models  of  the  lower  crust  and  upper  mantle. 
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OBJECTIVES 


We  describe  a  technique  based  on  waveform  fitting  by  synthetic  seismograms,  and  demonstrate  its  application  to 
determine  the  crust  and  upper  mantle  structure  beneath  the  continent  of  Africa.  Our  technique  utilizes  the  reflectivity 
method  (Kennett,  1 983)  to  compute  the  synthetic  seismograms  for  an  earthquake  recorded  at  a  particular  seismic 
station  and  implements  a  global  optimization  algorithm  using  VFSA  (Ingber,  1989;  Sen  and  Stoffa,  1995)  to  invert 
for  the  crust  and  upper  mantle  structure  beneath  the  station.  Our  technique  is  complimentary  to  the  existing  receiver 
function  analysis  method  (e  g.,  Owens  et  al.,  1987;  Ammon,  1991)  and  retains  its  advantages,  yet  minimizes  its 
potential  limitations  when  the  derived  structural  models  are  used  to  locate  and  determine  the  focal  depths  of  small, 
regional  events. 

RESEARCH  ACCOMPLISHED 


Modeling  Method 

Following  the  development  of  the  forward  problem  described  in  Pulliam  and  Sen  (2005),  we  perform  an 
optimization  procedure  to  determine  the  best  model  for  a  given  source- receiver  pair.  In  this  analysis  we  employ  a 
‘‘global  optimization  algorithm”  which  is  only  weakly  dependent  on  the  choice  of  the  initial  model.  In  particular,  we 
use  a  method  called  VFSA,  which  is  a  variant  of  simulated  annealing  (SA)  aimed  at  making  the  computations  more 
efficient  (Ingber,  1989;  Sen  and  Stoffa,  1995).  To  further  illustrate  the  VFSA  technique,  a  simplified  flow-chart  is 
shown  in  Figure  1 .  The  method  starts  with  an  initial  model  (m°)  with  an  associated  error  or  energy,  E(m°).  It  then 
draws  a  new  model,  mnew,  among  a  distribution  of  models  from  a  temperature  (T)  dependent  Cauchy-like 
distribution,  r(T),  centered  on  the  current  model  (Figure  1).  The  associated  error  or  energy,  E(mnew),  is  then 
computed  and  compared  with  E(m°)  (Figure  1).  If  the  change  in  energy  (6E)  is  less  than  or  equal  to  zero,  the  new 
model  is  accepted  and  replaces  the  initial  model.  However,  if  the  above  condition  is  not  satisfied,  mnew  is  accepted 
with  a  probability  of  [e6F >T]  (Figure  1).  This  rule  of  probabilistic  acceptance  in  SA  allows  it  to  escape  a  local 
minimum.  The  processes  of  model  generation  and  acceptance  are  repeated  a  large  number  of  times  with  the 
annealing  temperature  gradually  decreasing  according  to  a  predefined  cooling  schedule  (Figure  1).  VFSA  is  more 
efficient  than  the  traditional  SA  because  it  allows  for  larger  sampling  of  the  model  space  during  the  early  stages  of 
the  waveform  fitting,  and  much  narrower  sampling  in  the  model  space  as  the  procedure  converges  and  the 
temperature  decreases,  while  still  allowing  the  search  to  escape  from  the  local  mmima.  Additionally,  the  ability  to 
perform  different  perturbations  for  different  model  parameters  allows  for  individual  control  of  each  parameter  and 
the  incorporation  of  a  priori  information  (Sen  and  Stoffa,  1995). 

Solutions  to  geophysical  inverse  problems  are  often  nonunique.  It  is  therefore  necessary'  to  explore  the  model  space 
and  thus  identify  the  range  of  models  that  fit  the  data,  and  perhaps  to  identify  characteristics  of  the  models  that  are 
required  by  the  data,  rather  than  which  simply  are  allowed  by  the  data.  VFSA  conducts  such  a  search  efficiently,  and 
the  products  of  multiple  such  searches  enable  us  to  evaluate  the  uncertainty  in  a  single,  best-fitting  solution.  This 
evaluation  is  particularly  necessary  in  seismic  waveform  modeling  because  more  than  one  model  can  often  explain 
the  observed  data  equally  well,  and  trade-offs  between  different  model  parameters  are  common  (Pulliam  and  Sen, 
2005).  The  waveform  inversion  method  we  use  in  this  study  incorporates  important  statistical  tools  that  allow  the 
user  to  evaluate  the  uniqueness,  and  physical  feasibility  of  the  resulting  model.  The  most  useful  of  these  tools  in 
evaluating  the  results’  reliability  are  the  posterior  probability  density  (PPD)  function,  and  the  parameter  correlation 
matrix.  To  estimate  these  statistical  parameters  we  cast  the  inverse  problem  in  a  Bayesian  framework  (e  g..  Sen  and 
Stoffa,  1 995),  and  employ  “importance  sampling”  based  on  a  Gibbs’  sampler  (GS)  (Sen  and  Stoffa,  1995;  Pulliam 
and  Sen,  2005).  The  goal  of  “importance  sampling”  is  to  concentrate  sample  points  in  the  regions  that  are  the  most 
“significant,”  in  some  sense  (perhaps,  for  example,  where  the  error  function  is  rapidly  varying,  or  many  acceptable 
solutions  lie).  Because  this  concentration  is  achieved  using  a  Gibbs’  probability  distribution,  it  has  been  named  the 
"Gibbs’  sampler”  (Sen  and  Stoffa,  1995).  The  PPD  function  [a(m|d0^)]  is  defined  as  a  product  of  a  likelihood 
function  [e'£(m)],  and  prior  probability'  density  function, p(m).  The  prior  probability  density  function  p{ m),  describes 
the  available  infonnation  on  the  model  without  the  knowledge  of  the  data  and  defines  the  probability  of  the  model  m 
independent  of  the  data.  The  likelihood  function  defines  the  data  misfit  and  its  choice  depends  on  the  distribution  of 
error  in  the  data  (Sen  and  Stofta,  1 995).  Sen  and  Stofta  (1996)  examined  several  different  approaches  to  sampling 
models  from  the  PPD  and  concluded  that  a  multiple- VFSA  based  approach,  though  theoretically  approximate,  is  the 
most  efficient.  Here,  we  have  employed  the  fast  multiple- VFSA  algorithm  to  estimate  uncertainties.  We  compute 
approximate  marginal  PPD  and  posterior  correlation  matrices  to  characterize  uncertainties  in  the  derived  results.  The 
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posterior  correlation  matrix  measures  the  relative  trade-off  between  individual  model  parameters  and  is  computed  by 
normalizing  the  covariance  between  two  model  parameters  (Sen  and  Stoffa,  1996). 

Application 

We  apply  the  modeling  method  outlined  above  to  data  from  large-magnitude,  deep-focus  earthquakes  recorded 
teleseism  ically  at  twelve  permanent  broadband  seismic  stations  spanning  the  continent  of  Africa  (Figure  2).  A  total 
of  seventeen  earthquakes  are  used  in  this  study  selected  from  the  global  catalog  of  Harvard  Centroid  Moment 
Tensors  (CMT)  (1976-2004)  (Table  1).  Seven  of  these  are  located  in  the  Hindu  Kush,  three  in  southern  Bolivia,  two 
each  in  Argentina,  and  the  Java  Sea,  and  one  each  in  Afghanistan,  South  Sandwich  Islands,  and  southern  Sumatra. 
The  focal  depths  of  these  earthquakes  range  between  200  km-600  km,  and  their  magnitudes  lie  between  5. 5-7.0. 
Since  the  goal  of  this  study  is  to  also  model  the  SPL  phase  that  is  generated  close  to  the  seismic  station,  we  choose 
such  a  focal  depth  range  to  eliminate  the  (SPL)  phase  generated  at  the  earthquake  source.  Epicentral  distances  from 
the  seismic  stations  of  the  selected  earthquakes  are  between  30°  and  80°  so  as  to  avoid  possible  incorporation  of 
phases  that  interacted  with  the  earth’s  core.  Based  on  our  selection  criteria,  earthquake  data  available  from  the 
seismic  stations  in  Africa  are  sparse.  Hence,  the  numbers  of  earthquakes  recorded  by  each  station  range  from  1  to  9 
and  azimuthal  coverage  is  poor.  In  spite  of  the  sparse  set  of  source-receiver  pairs  available  for  modeling,  we  are  able 
to  obtain  reliable  azimuthally  dependent  models  in  some  cases. 


Table  1.  List  of  earthquakes  used  in  this  study 


Event 

Number 

Date 

(>yyy  mm  dd) 

Time 

(hr:  min:  sec) 

Latitude 

Longitude 

Focal  Depth 
(kin) 

Magnitude 

Region 

1 

1990  07  13 

14:27:24  S 

36  46  N 

70  80  E 

218 

5.6 

Hindu  Kush 

2 

1991  06  23 

21:31:58.33 

26.82  S 

63.40  W 

581 

6  4 

Santiago  del 
Estero.  Argentina 

3 

1993  08  09 

11:44:08.306 

36  42  N 

70.72  E 

210 

5.8 

Hindu  Kush 

4 

1993  08  09 

12  48  29  404 

36  36  N 

70  85  E 

230 

6  3 

Hindu  Kush 

5 

1994  06  30 

09:31:04  212 

36.25  N 

71  08  E 

233 

6  1 

Afghanistan 

6 

1994  10  25 

00:58:47.987 

36.30  N 

70.91  E 

244 

5.9 

Hindu  Kush 

7 

1994  11  15 

20:27:46.338 

5  61  S 

1 10.2  E 

559 

6.2 

Java  Sea 

8 

1994  12  07 

03:47:10.307 

23.46  S 

66  74  W 

243 

5  6 

Jujuy  Province. 
Argentina 

9 

1995  05  13 

21:07:41.55 2 

5.22  S 

108.92  E 

554 

5.7 

Java  Sea 

10 

1997  01  23 

02:24:53.890 

22  S 

65.72  W 

276 

64 

Southern  Bolivia 

11 

1997  10  05 

18:14:56.228 

59.74  S 

29.20  W 

273.9 

6.0 

South  Sandwich 
Islands 

12 

1997  12  17 

05:55:49.933 

36.39  N 

70.77  E 

207 

5.5 

Hindu  Kush 

13 

1998  03  21 

18:26:09.8 

36  43  N 

70  13  E 

227.8 

5.8 

Hindu  Kush 

14 

1999  09  15 

03:11:17.074 

20.93  S 

67.28  W 

218 

60 

Southern  Bolivia 

15 

2002  03  03 

12:12:00.183 

36  43  N 

70.44  E 

209 

6  3 

Hindu  Kush 

16 

2003  07  27 

11  49:52.347 

20.13  S 

65.18  W 

345.3 

5.9 

Southern  Bolivia 

17 

2004  07  25 

14:44:40  148 

2  43  S 

103  93  E 

582.1 

68 

Southern  Sumatra 
Indonesia 

Results  of  this  Application 

Below,  we  report  results  of  waveform  fitting  for  selected  teleseismic  earthquakes  recorded  in  Africa.  A  comment  on 
amplitude  matches:  The  most  successful  match  between  synthetics  and  data  would  be  one  in  which  the  synthetic 
waveform  matched  the  data  exactly — wiggle  for  wiggle.  This  is  unrealistic  for  several  reasons,  including  the  fact 
that  models  used  to  compute  synthetics  are  layered,  isotropic,  limited  to  ten  to  sixteen  layers,  and  have  fixed 
attenuation  (Q)  values.  Further,  the  source  time  function  is  assumed  to  be  Gaussian  and  its  focal  mechanism  is 
assumed  to  be  correctly  represented  by  Harvard’s  CMT  solution.  To  minimize  complexities  in  the  source  time 
function  we  avoid  very  large  earthquakes.  Given  the  uncertainties  in  model  Q  and  focal  mechanisms,  which  will 
largely  control  relative  amplitudes  of  various  phases,  we  focus  our  fitting  criteria  on  matching  each  phase’s  arrival 
time,  polarity,  and  pulse  character.  Fitting  the  amplitude  of  each  phase,  while  desirable,  is  deemed  to  be  of  lesser 
importance. 
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West  Africa 

The  region  encompassing  north  and  west  Africa  includes  the  seismic  stations  of  TAM,  DBIC,  MBO,  and  MDT 
(Figure  2).  Among  these  stations,  TAM  recorded  data  of  better  quality;  examples  of  the  waveform  fits  for  events  1 , 
3,  and  6  (Table  1)  recorded  at  TAM  are  shown  in  Figure  3a.  We  observe  S,  SP,  and  SsPmP  phases  consistently  on 
all  these  event  seismograms,  and  they  correlate  well  with  the  synthetics  generated  by  the  optimization  technique 
(Figure  3a).  On  event  3,  we  also  observe  a  prominent  SPL  phase  and  the  synthetics  match  it  well  (Figure  3a). 

Particle  motion  diagrams  for  the  corresponding  time  window  on  both  data  and  synthetics  confirm  this  observation 
(Figure  3b).  We  observe  prograde  elliptical  particle  motion,  which  is  diagnostic  of  the  SPL  phase,  on  both  diagrams 
(Figure  3b).  Except  for  event  3,  none  of  the  others  have  any  signature  of  the  SPL  phase,  as  confirmed  by  particle 
motion  diagrams  for  corresponding  time  windows.  This  observation — that  one  source-receiver  pair  would  show  SPL 
while  other,  similar  paths  would  not — is  unexpected  and  we  cannot  explain  it. 


Figure  1.  Flow-chart  elaborating  the  very  fast  simulated  annealing  (VFSA)  algorithm  used  in  this  study  for 
the  waveform  inversion  by  global  optimization  (Modified  from  Sen  and  Stoffa,  1995).  E(m°)-error 
function  for  the  initial  model  m°,  E(mnew)-error  function  for  the  new  model  m1*",  T-temperature, 
r(T)-temperature  dependent  Cauchy-like  distribution. 
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Figure  2.  Map  showing  locations  of  the  earthquakes  used  in  this  study  (blue  stars)  and  the  permanent 
broadband  seismic  stations  in  the  African  continent  (red  triangles)  that  have  recorded  these 
earthquakes.  The  respective  station  codes  are  shown  adjacent  to  each  seismic  station. 


Based  on  the  waveform-fitting  results  for  all  six  events  recorded  at  TAM,  we  generate  P-  and  S-wave  velocity 
models  up  to  a  depth  of  100  km  (Figure  3c).  We  observe  some  variability  in  the  models  (Figure  3c),  and  hence 
compute  the  uncertainties  for  each  model  using  the  statistical  tools  described  earlier  to  choose  the  “best”  model.  In 
Figure  3d,  we  show  examples  of  parameter  correlation  matrices  computed  from  the  modeling  results  of  events  1  and 
3  (Table  1 ).  Each  small  square  along  an  axis  of  the  parameter  correlation  matrix,  either  horizontally  or  vertically, 
represents  a  model  parameter  (Figure  3d).  Since  every  model  layer  consists  of  four  independent  model  parameters 
(Vp,  Vs,  thickness,  and  density),  four  small  squares  combined  together  represent  a  model  layer  on  both  axes 
(Figure  3d).  Correlation  values  range  between  -1  and  1  and  are  symmetric  about  the  diagonal  of  the  matrix,  hence, 
for  clarity,  we  show  only  values  below  the  diagonal  (Figure  3d).  Values  along  the  diagonal  are  ones,  simply 
indicating  that  each  parameter  is  perfectly  correlated  with  itself.  Oft -diagonal  colored  squares  indicate  significant 
cross-correlation  (trade-offs)  between  corresponding  model  parameters.  In  the  parameter  correlation  matrices  for 
both  events  (Figure  3d),  layers  comprising  the  upper  crust  have  greater  independence,  as  indicated  by  the  sparse 
distribution  of  off-diagonal  cross-correlations  whose  absolute  values  are  greater  than  ±0.5  (colored  squares). 

Also,  for  both  events,  the  level  of  tradeoffs  among  model  parameters  in  these  shallow  layers  is  similar 
(Figure  3d).  For  event  1 ,  however,  the  layers  comprising  the  lower  crust  and  upper  mantle  have  larger  off-diagonal 
cross-correlations,  indicating  significant  tradeoffs  (Figure  3d).  On  the  contrary,  for  event  3,  even  the  lower  crustal 
and  upper  mantle  layers  appear  better  constrained  (Figure  3d).  Intriguingly,  the  SPL  phase  is  also  observed  in  the 
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seismogram  of  event  3  but  not  in  that  of  event  1  (Figure  3a).  This  observation  attests  to  the  fact  that,  if  SPL  is 
present  in  the  seismogram  and  is  well  modeled,  we  are  able  to  better  constrain  the  structure  of  the  lower  crust  and 
upper  mantle.  This  result,  which  is  expected,  due  to  the  sensitivity  of  SPL  to  those  parts  of  the  model  (Figure  1  b), 
drives  our  decision  to  generate  velocity  models  down  to  the  Moho  for  seismic  stations  at  which  SPL  is  not  observed 

We  observe  and  successfully  fit  S  and  SP  phases  in  seismograms  recorded  at  DBIC  and  at  MBO.  The  SsPmP  phase 
appears  on  the  vertical-component  seismograms  of  events  13  at  DBIC  and  event  16  at  MBO  but  is  absent  on  their 
radial-component  seismograms.  We  do  not  observe  the  SPL  phase  in  seismograms  recorded  at  DBIC  but  do  so  at 
MBO  in  the  seismograms  of  event  16.  The  SPL  phase  is  prominent  on  both  vertical  and  radial  component 
seismograms  of  event  16  recorded  at  MBO,  which  we  confirm  with  particle  motion  diagrams  that  show  prograde 
elliptical  motion  for  the  corresponding  time  window.  At  MDT,  within  the  time-window  expected  to  contain  the 
phases  analyzed  in  this  study,  we  are  unable  to  clearly  identify  them  and  they  appear  to  be  contaminated  by 
interfering  arrivals,  and  hence  are  also  not  well- correlated  with  the  synthetics  generated  by  the  waveform  inversion 
process.  Station  MDT  is  located  near  the  Atlas  Mountains  in  Morocco,  and  so  may  be  underlain  by  complicated 
three-dimensional  structure,  which  the  waveform  modeling  program  used  in  this  study  is  unable  to  model  accurately 
The  lack  of  correlation  between  synthetic  seismograms  and  data  at  MDT  is  likely  a  consequence  of  this  limitation 

East  Africa 

Stations  ATD,  FURI,  KMBO,  and  MBAR  are  situated  within  and  on  the  flanks  of  the  active  East  African  rift  system 
and  waves  recorded  at  these  stations  sample  the  complicated  three-dimensional,  anisotropic  structure  beneath  the  rift 
(Ayele  et  al.,  2004;  Dugda  and  Nyblade,  2006).  Three-dimensional  structure  is  manifested  in  the  seismograms  as 
numerous,  possibly  scattered,  refracted,  or  split,  phase  arrivals  with  strong  interference  amongst  themselves. 
Anisotropy  inferred  from  shear-wave  splitting  studies  has  also  been  reported  for  stations  ATD,  FURI,  and  KMBO 
by  Ayele  et  al.  (2004).  The  waveform  inversion  method  we  use  in  this  study  is  capable  of  estimating  only 
one-dimensional  (although  azimuthally  dependent)  structure  only,  hence  waveforms  for  events  recorded  at  these 
stations  are  not  precisely  correlated  in  some  cases.  On  both  the  vertical  and  radial  component  seismograms  of  most 
events  at  ATD,  FURI,  KMBO,  and  MBAR,  we  observe  S  and  SsPmP  phases.  On  the  other  hand,  except  at  MBAR, 
we  observe  the  SP  phase  only  on  the  vertical  component  for  these  events.  Only  at  FURI,  for  event  1 2  (Table  1 ),  and 
MBAR,  for  event  1 5  (Table  1),  do  we  see  SPL  phases  in  the  seismograms. 

Southern  Africa 

Seismic  stations  TSUM,  LSZ,  LBTB,  and  SUR  are  located  in  southern  Africa  (Figure  2).  However,  for  the  lone 
event  recorded  at  SUR,  we  are  able  to  identify  only  the  direct  S  phase,  and  thus  we  do  not  attempt  to  generate  P-  and 
S-wave  velocity  models  for  SUR.  We  observe  and  successfully  fit  S,  SP,  and  SsPmP  phases  on  events  13  and  14 
(Table  1)  recorded  at  TSUM,  event  9  (Table  1)  recorded  at  LSZ,  and  event  1 1  (Table  1)  at  LBTB.  But,  except  for 
event  1 3  at  TSUM,  we  are  able  to  identify  the  SP  phase  on  both  the  vertical  and  radial  component  seismograms  at 
TSUM,  LSZ,  and  LBTB.  Similarly,  the  SsPmP  phase  is  only  identifable  on  the  vertical  component  seismogram  for 
events  14  and  9  recorded  at  TSUM  and  LSZ  respectively.  It  is,  however,  absent  on  the  seismograms  for  event  1 3  at 
TSUM.  We  do  not  observe  the  SPL  phase  in  the  seismograms  of  any  event  recorded  at  these  stations. 

Summary  of  Models 

Figures  4a-d  show  P-  and  S-wave  velocity  models  for  the  north  and  west  African  stations  of  TAM,  DBIC,  MBO, 
and  MDT.  Estimates  of  crustal  thickness  beneath  these  stations  range  between  36  km  and  42  km  (Figures  4a-d), 
comparable  to  regional  estimates  of  34  km  to  40  km  by  Pasyanos  et  al.  (2004).  We  also  note  that  the  crust  is  slightly 
thicker  in  west  Africa,  e  g.  beneath  stations  DBIC  and  MBO  (--41-42  km)  (Figures  4b  and  4c),  compared  to  the 
seismic  stations  TAM  and  MD I  in  north  Africa  (~36  km— 38  km)  (Figures  4a  and  4d).  A  similar  observation  was 
also  made  earlier  by  Pasyanos  and  Walter  (2002)  using  surface  wave  dispersion  tomography,  and  by  Marone  et  al. 
(2003)  using  joint  inversion  of  local,  regional,  and  teleseismic  data.  Except  for  TAM,  the  crust  below  all  the  stations 
appears  to  be  fairly  simple  in  structure  (Figures  4a-d)  (Sandvol  et  al.,  1 998),  suggesting  that  it  is  minimally  affected 
by  large-scale  tectonic  processes.  However,  a  middle  to  lower  crustal  low-velocity  zone  obtained  beneath  all  the 
seismic  stations  in  the  region  (Figures  4a-d),  indicate  possible  local  tectonic  influences.  Our  estimate  of  crustal 
thickness  beneath  TAM  (~36  km)  is  similar  to  that  obtained  from  receiver  function  studies  by  Sandvol  et  al.  (1998) 
(38  ±  2  km)  (Figure  4a),  from  Rayleigh  wave  group  velocity  dispersion  studies  by  Hazier  et  al.  (2001 )  (43  ±  5  km), 
and  from  surface- wave  dispersion  tomography  by  Pasyanos  and  Walter  (2002)  (~40  km).  TAM  is  close  to  the 
location  of  the  Hoggar  hot  spot  but,  as  noted  by  Sandvol  et  al.  ( 1 998),  the  crustal  thickness  indicates  that  a  mantle 
plume  has  not  significantly  altered  the  crust  here.  However,  in  contrast  to  the  model  of  Sandvol  et  al.  (1998),  the 
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crustal  P-  and  S-wave  velocities  we  obtain  in  this  study  at  TAM  are  both  slightly  lower  (Figure  4a).  These  velocities 
at  TAM  range  between  6.25  km/s-6.8  km/s,  and  3. 1  km/s-3.9  km/s,  respectively  (Figure  4a).  Upper  mantle  P-wave 
velocities  exhibit  a  gradational  increase  with  depth  below  the  Moho,  whereas  the  S-wave  velocities  are  nearly 
constant  (4  km/s-4.2  km/s)  within  that  range  of  depths  (Figure  4a).  Furthermore,  we  also  obtain  an  anomalous  low- 
velocity  zone  of~5  km  thickness  at  the  base  of  the  upper  crust  (Figure  4a)  that  appears  to  be  well  constrained  based 
on  the  uncertainty  estimates  described  earlier  (Figure  3d).  However,  we  are  unable  to  confirm  the  existence  of  this 
layer  from  any  other  independent  studies. 


(a) 


(b) 


Figure  3.  (a)  Vertical  and  radial  component  seismograms  for  events  1,  3,  and  6  (Table  1)  recorded  at  TAM 
showing  the  observ  ed  (solid  line)  and  synthetic  (dashed  line)  waveforms.  The  correlated  w  aveforms 
are  indicated  on  the  panels,  (b)  Particle  motion  diagrams  for  the  portion  of  the  w  aveforms  in  the  time 
window  (1077—1085  s)  on  the  data  and  synthetics  for  event  3  showing  prograde  ellipitical  motion 
diagnostic  of  the  SPL  phase.  The  dotted  portions  of  the  diagrams  indicate  beginning  of  the  motion. 
Figure  3  is  continued  on  next  page. 
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Figure  3  continued,  (c)  P-  and  S-wave  velocity  models  up  to  100  km  for  station  TAM  from  the  inversion 
results  for  individual  events  recorded  at  TAM.  (d)  Model  parameter  correlation  matrices  for 
events  1  (left  panel)  and  3  (right  panel).  Each  small  square  represents  a  model  parameter  (Vp,  Vs, 
thickness  of  layer,  and  density)  on  both  the  axes.  The  correlations  range  between  -1  and  1.  Sparse 
off-diagonal  colored  squares  in  the  lower  crust-upper  mantle  in  event  3  (right  panel)  compared  to 
that  in  event  1  (left  panel)  indicate  better  resolution  and  confidence  (less  tradeoff)  in  this  region. 
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Figure  4.  P-  and  S-wave  velocity  models  (solid  lines)  obtained  for  (a)  TAM  (b)  DBIC  (c)  MBO  (d)  MDT 
(e)  ATD  (f)  FURI  (g)  KMBO  (h)  MBAR  (i)  TSUM  (j)  LSZ  (k)  LBTB.  The  P-  and  S-wave  velocity 
models  (broken  lines)  in  (a),  (b),  (d),  and  (e)  are  from  receiver  function  studies  by  Sandvol  et  al. 
(1998),  in  (0  are  from  receiver  function  studies  by  Ayele  et  al.  (2004),  and  in  (j)  and  (k)  are  from 
receiver  function  studies  by  Midzi  and  Ottemoller  (2001). 
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At  the  west  African  coastal  station  DBIC,  we  obtain  a  Moho  depth  of  -41  km  which  is  similar  to  that  obtained  by 
Sandvol  et  al.  (1998)  (-40  ±  2.3  km)  (Figure  4b).  P-wave  velocities  range  between  6.7  km/s-7  km/s  in  the  upper 
crust,  and  6.5  km/s-7. 3  km/s  in  the  lower  crust  (Figure  4b).  On  the  contrary,  S-wave  velocities  show  a  gradational 
increase  in  the  crust  with  depth  from  3.7  km/s^l.7  km/s  (Figure  4b).  Anomalous  Vp/Vs  ratios  are  thus  caused  by  a 
low  P-wave  velocity  zone  of  -15  km  thickness  in  the  lower  crust.  However,  due  to  the  trade-offs  between  the  model 
parameters  in  this  depth  range,  we  conclude  that  this  anomaly  is  not  well-constrained.  Beneath  MBO,  any  prior 
P-  and/or  S-wave  velocity  models  are  absent.  Sandvol  et  al.  (1998)  had  analyzed  seismic  data  recorded  at  MBO  but 
the  anomalous  data  did  not  allow  them  to  estimate  a  velocity  model.  We  obtain  a  crustal  thickness  of -42  km 
(Figure  4c),  which  is  similar  to  that  found  underneath  other  stations  in  the  region.  A  regional  crustal  thickness  of 
43  ±  5  km  obtained  from  Rayleigh  wave  group  velocity  dispersion  studies  (Hazier  et  al.,  2001)  correlates  well  with 
the  results  of  this  study  at  MBO.  The  P-  and  S-wave  velocities  at  MBO  in  the  crust  range  between  5.6  km/s-7. 2 
km/s,  and  3.1  km/s-4. 1  km/s,  respectively  (Figure  4c).  We  also  observe  an  anomalous  lower  crustal,  approximately 
15  km-thick  zone  of  relatively  low  P-  and  S-wave  velocities  (6.8  km/s  and  3.8  km/s)  beneath  MBO  (Figure  4c).  Our 
P-  and  S-wave  velocity  models  beneath  MDT  predict  a  Moho  depth  of  -38  km  (Figure  4d).  Surprisingly,  even  with 
the  poor  waveform  fit  by  synthetics  to  the  event  recorded  at  MDT,  this  result  is  consistent  with  the  estimates 
obtained  by  Sandvol  et  al.  (1998)  (36  ±  1 .3  km).  As  also  noted  by  Sandvol  et  al.  (1998)  and  Pasyanos  and  Walter 
(2002),  the  slightly  shallower  Moho  at  MDT,  compared  to  that  at  DBIC  and  MBO,  indicates  that  in  spite  of  its 
proximity  to  the  Atlas  Mountains,  there  is  no  crustal  thickening  associated  with  them,  and  a  significant  root  is  absent 
beneath  the  mountains.  Sandvol  et  al.  (1998)  concluded  that  this  may  be  a  possible  outcome  of  the  fact  that  there 
existed  a  failed  rift  earlier  which  was  subsequently  inverted.  In  this  study,  beneath  MDT,  we  obtain  average  P-  and 
S-wave  velocities  in  the  crust  of  -6.4  km/s  and  3.7  km/s,  respectively,  except  in  a  low-velocity  zone  of -10  km 
thickness  in  the  lower  crust  where  these  are  6.2  km/s  and  3.4  km/s,  respectively  (Figure  4d).  A  similar  zone  has  also 
been  postulated  by  the  earlier  velocity  model  obtained  from  receiver  function  studies  (Sandvol  et  al,  1998). 
However,  because  of  poor  waveform  fits  at  MDT  in  this  study,  we  are  unable  to  postulate  the  existence  of  this  zone. 

For  East  Africa,  we  generate  P-  and  S-wave  velocity  models  beneath  seismic  stations  ATD,  FURI,  KMBO,  and 
MBAR,  which  are  shown  in  Figures  4e-h.  Located  within  and  on  the  flanks  of  the  East  African  Rift  System,  which 
is  relatively  well  studied,  these  stations  are  situated  in  active  tectonic  domains.  Except  beneath  ATD  (Figure  4e), 
where  the  crust  is  significantly  thin  compared  to  the  other  stations  in  the  region  (Figures  4f-h),  the  Moho  is 
generally  between  -38  km-41  km  deep.  The  estimate  of  crustal  thickness  beneath  ATD,  however,  is  the  subject  of 
an  active  debate.  Using  a  grid  search  method  to  model  receiver  functions  for  eleven  earthquakes  recorded  at  ATD, 
Sandvol  et  al.  (1998)  obtained  a  crustal  thickness  of -10  km  (Figure  4e).  But,  Dugda  and  Nyblade  (2006)  used  H-k 
analysis  of  receiver  functions  and  predicted  a  crustal  thickness  of  -23  ±  1 .5  km  beneath  ATD,  consistent  with  earlier 
results  from  inversion  of  gravity  data  for  the  general  area  by  Tiberi  et  al.  (2005).  In  this  study,  our  velocity’  model 
beneath  ATD  shows  comparable  velocity  discontinuities  at  -10  km  and  -21  km  depths  (Figure  4e),  suggesting  that 
either  of  these  depths  could  be  interpreted  as  the  Moho.  However,  the  layer  at  10  km  depth  appears  to  be  poorly 
constrained  compared  to  the  layer  at  —21  km  depth,  as  evidenced  from  the  PPD  and  the  parameter  correlation  matrix 
computations.  Therefore,  we  prefer  a  crustal  thickness  of  -21  km.  Irrespective  of  the  debate  on  the  crustal  thickness 
at  ATD,  the  crust  beneath  it  is  significantly  thinner  than  the  crust  beneath  other  seismic  stations  in  east  Africa 
(Figures  4f-h).  Located  within  the  Afar  depression,  close  to  the  coast  of  the  Red  Sea  on  the  eastern  edge  of  the 
African  continent,  such  a  thin  crust  is  expected  at  ATD  because  of  highly  stretched  continental  crust  (Sandvol  et  al, 
1998).  The  P-  and  S-wave  velocities  in  the  crust  at  ATD  range  between  4.7  km/s-7.2  km/s,  and  2.5  km/s-4. 3  km/s, 
respectively  (Figure  4e).  Crustal  P-wave  velocities  below  about  5  km  depth  are  relatively  high  and,  as  noted  by 
Dugda  and  Nyblade  (2006),  could  mdicate  a  highly  mafic  composition  caused  by  igneous  rock  emplacement  during 
the  syn-rift  stage.  Figure  4f  shows  the  preferred  P-  and  S-wave  velocity  models  beneath  FURI  which  is  situated  in 
the  northern  part  of  the  western  Ethiopian  plateau.  We  obtain  an  estimate  of  crustal  thickness  beneath  FURI  of 
-39  km,  which  is  similar  to  that  obtained  using  receiver  function  analyses  by  Ayele  et  al.  (2004)  (-40  km),  and 
Dugda  et  al.  (2005)  (-44  km).  The  -40  km  thick  crust  beneath  FURI,  which  is  located  close  to  the  border  of  the 
western  Ethiopian  plateau  and  the  Afar  depression,  is  also  consistent  with  previous  refraction  studies  of  Berckhemer 
et  al.  (1975)  as  reported  by  Ayele  et  al.  (2004).  Crustal  P-  and  S-wave  velocities  below  -5  km  depth  at  FURI  range 
between  5.3  km/s-6.8  km/s,  and  3.2  km/s-3.6  km/s,  respectively  (Figure  4f).  In  addition,  beneath  FURI,  our  results 
also  predict  an  -50  km  thick  layer  immediately  below  the  Moho  in  the  upper  mantle  that  has  P-  and  S-wave 
velocities  of -7. 1  km/s  and  -4.3  km/s,  respectively,  which  are  anomalously  slow  (Figure  4f).  A  similar  layer  with  P- 
and  S-wave  velocities  of -7.4  km/s  and  -4.2  km/s,  respectively,  was  also  obtained  by  Ayele  et  al.  (2004)  (Figure 
4f).  As  noted  by  Ayele  et  al.  (2004),  this  anomalously  slow  layer  possibly  indicates  altered  lithospheric  material,  and 
supports  the  result  from  Rayleigh  wave  dispersion  by  Knox  et  al.  (1999),  of  an  approximately  100  km  thick 
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lithosphere  beneath  FURL  Station  KMBO  is  located  in  Kenya,  close  to  the  southern  end  of  the  eastern  branch  of  the 
East  African  Rift  System,  but  outside  its  edge.  Beneath  KMBO  our  estimate  of  the  crustal  thickness  is  -38  km 
(Figure  4g).  This  estimate  is  similar  to  that  obtained  using  receiver  function  analysis  by  Dugda  et  al.  (2005)  (-41 
km).  Crustal  P-  and  S-wave  velocities  show  a  gradational  increase  with  depth  and  range  between  5.4  km/s-8  km/s, 
and  3.5  km/s— 4.5  km/s,  respectively  (Figure  4g).  The  velocity  structure  beneath  KMBO  appears  to  be  fairly  simple. 
Although  it  has  relatively  high  P-wave  velocities  in  the  lower  crust  (Figure  4g),.it  is  otherwise  typical  of  cratonic 
regions.  Seismic  station  MBAR  is  located  between  the  western  boundary  of  the  Tanzania  craton  and  the  western 
branch  ol  the  East  African  Rift  System.  Due  to  poor  correlation  of  some  ot  the  observed  phases  with  synthetics,  and 
the  availability  of  only  one  event,  the  velocity  model  beneath  MBAR  is  poorly  constrained.  Nevertheless,  our 
estimate  of  crustal  thickness  beneath  MBAR,  to  our  knowledge  the  first  of  its  kind,  is  -41  km  (Figure  4h),  and  is 
consistent  with  that  obtained  from  other  stations  in  the  region.  Crustal  P-  and  S-wave  velocities  at  MBAR  range 
between  5.3  km/s-7.7  km/s,  and  3.2  km/s-3.8  km/s,  respectively  (Figure  4h).  Our  model  also  predicts  a  low  P-wave 
velocity  (-7.2  km/s)  layer  beneath  the  crust  (Figure  4h).  Such  a  layer  promotes  the  generation  of  the  SPL  phase  near 
the  station,  which  we  observe  in  both  data  and  synthetics  for  the  event  recorded  at  MBAR.  Therefore,  in  spite  of  the 
poorly  constrained  model  obtained  in  this  study,  we  cannot  rule  out  the  possibility  of  its  existence. 

In  southern  Africa  we  generated  P-  and  S-wave  velocity  models  beneath  the  seismic  stations  TSUM,  LSZ,  and 
LBTB  (Figures  4i-k).  Although  we  analyzed  seismic  data  recorded  at  SUR,  due  to  the  lack  of  identifiable  phases, 
we  do  not  generate  P-  and  S-wave  velocity  models  for  the  station.  Similar  to  most  of  the  results  in  north  and  west 
Africa,  our  results  for  southern  Africa  are  representative  of  stable  shield  regions.  However,  in  general,  we  obtain 
slightly  higher  crustal  thicknesses  ranging  between  -42  km  and  46  km  (Figures  4i-k).  We  also  predict  crustal  low 
velocity  zones  as  discussed  later  in  our  models  beneath  two  of  the  three  stations  in  southern  Africa.  Beneath  TSUM, 
to  our  knowledge,  no  prior  velocity  model  exists.  Thus,  the  velocity  model  obtained  from  our  study  at  TSUM  is  the 
first  of  its  kind.  We  obtained  a  crustal  thickness  beneath  TSUM  of -42  km  (Figure  4i).  The  velocities  of  P-  and  S- 
waves  in  the  crust  range  between  6.3  km/s-7.3  km/s,  and  3.2  km/s-4  km/s,  respectively  (Figure  4i).  These  results 
are  similar  to  those  obtained  for  the  seismic  stations  in  north  and  west  Africa  and  are  therefore  representative  of 
stable  shield  regions.  We  do  not  obtain  any  anomalous  P-  and  S-wave  velocity  zones  beneath  TSUM  (Figure  4i).  At 
LSZ,  our  study  indicates  that  the  Moho  is  located  at  a  depth  of  -43  km  (Figure  4j),  which  is  consistent  with  that 
obtained  by  Midzi  and  OttemOller  (2001)  (-40-43  km).  The  crustal  P-wave  velocity  is  nearly  constant  (-6.2  km/s) 
between  —8  km  to  32  km  depth  (Figure  4j).  Below  —32  km  depth,  P-wave  velocities  increase  rather  sharply  from 
-6.2  km/s  to  7 .8  km/s  at  the  Moho  (Figure  4j).  The  crustal  S-wave  velocities  range  between  -3.6  km/s  and  3.8  km/s 
(Figure  4j).  We  also  obtain  a  lower  crustal  low-velocity  zone  of -5-8  km  thickness  in  our  velocity  model  for  LSZ 
(Figure  4j),  consistent  with  models  for  seismic  stations  elsewhere  in  Africa  in  the  cratonic  regions.  Such  a 
phenomenon  was  also  noted  by  Midzi  and  OttemOller  (2001).  Figure  4k  shows  the  P-  and  S-wave  velocity  models 
beneath  LBTB  obtained  in  our  study.  There  appears  to  be  a  broad  crust-mantle  transition  zone  beneath  LBTB  and 
the  upper  bound  of  the  estimate  of  crustal  thickness  beneath  LBTB  is  -46  km  (Figure  4k).  Midzi  and  OttemOller 
(2001 )  also  noted  the  same  and  predicted  a  crust-mantle  transition  zone  between  37-45  km.  The  crustal  P-  and 
S-wave  velocities  beneath  LBTB  range  between  5.8  km/s— 7.5  km/s,  and  3.5  km/s— 4.2  km/s,  respectively,  except  for 
a  distinct  low  P-velocity  (5.4  km/s)  zone  of  -8  km  thickness  in  the  upper  crust  between  10  km-20  km  depth  (Figure 
4k).  Such  a  low-velocity  zone  was  also  obtained  by  Midzi  and  OttemOller  (2001)  at  similar  depths,  however,  the 
P-wave  velocities  predicted  from  our  study  for  this  zone  are  significantly  lower  than  those  predicted  by  Midzi  and 
OttemOller  (2001).  Given  its  appearance  beneath  other  cratonic  seismic  stations  in  Africa,  the  crustal  low-velocity 
zone  appears  to  be  a  general  characteristic  of  the  region. 


CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  paper,  we  discuss  a  waveform  fitting  technique  that  relies  on  a  parallelized  reflectivity  method  to  compute 
synthetic  seismograms  and  implements  a  global  optimization  algorithm  using  VFSA.  We  also  demonstrate  the 
application  of  the  method  to  determine  one-dimensional,  azimuthally  dependent,  crust  and  upper  mantle  P-  and 
S-wave  velocity  structure  beneath  broadband  seismic  stations  across  the  continent  of  Africa.  Our  technique  avoids 
dependence  of  the  final  results  on  the  initial  model,  and  we  are  able  to  compute  synthetic  seismograms  that  contain 
all  the  possible  phases  for  a  prescribed  source-receiver  path,  and  obtain  direct  estimates  of  the  P-  and  S-wave 
velocities  beneath  seismic  stations.  Statistical  tools  incorporated  in  the  technique  allow  us  to  assess  uncertainties 
associated  with  our  models  and  estimate  tradeoffs  between  model  parameters  in  different  layers.  The  use  of  the 
SPL  phase  as  shown  in  the  study,  enhances  our  constraints  for  lower  crust  and  upper  mantle  structure  beneath  the 
seismic  stations. 
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Applied  to  large-magnitude,  deep-focus  earthquakes  recorded  teleseismically  in  Africa,  our  method  successfully 
produced  crust  and  upper  mantle  (wherever  SPL  was  observed)  P-  and  S-wave  velocity  models,  that  are  consistent 
with  earlier  models,  in  the  sense  that  they  fall  within  the  associated  uncertainties  we  found  with  the  products  of 
multiple  VFSA  runs.  For  some  seismic  stations,  our  study  provided  such  velocity  models  that  are  the  first  of  their 
kind.  Our  models  were  also  consistent  with  the  regional  tectonics  of  Africa.  While  the  technique  described  here 
provided  layered,  one-dimensional  models,  a  dataset  that  includes  a  broader  azimuthal  distribution  of  earthquakes 
for  each  station  would  allow  this  source-receiver-based  technique  to  produce  better  azimuthal ly-dependent  models, 
and  thus  a  more  detailed  view  of  the  earth  s  structure.  We  plan  to  apply  this  method  to  determine  the  crust  and  upper 
mantle  structure  beneath  China  and  the  Middle  East. 
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